{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using ChoiceRank to understand network traffic\n",
    "\n",
    "This notebook provides a quick example on how to use ChoiceRank to estimate transitions along the edges of a network based only on the marginal traffic at the nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import choix\n",
    "import networkx as nx\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Generating sample data\n",
    "\n",
    "First, we will generate sample data.\n",
    "This includes\n",
    "\n",
    "1. generating a network,\n",
    "2. generating a parameter for each node of the network,\n",
    "3. generating samples of choices in the network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_items = 8\n",
    "p_edge = 0.3\n",
    "n_samples = 3000\n",
    "\n",
    "# 1. Generate a network.\n",
    "graph = nx.erdos_renyi_graph(n_items, p_edge, directed=True)\n",
    "\n",
    "# 2. Generate a parameter for each node.\n",
    "params = choix.generate_params(n_items, interval=2.0)\n",
    "\n",
    "# 3. Generate samples of choices in the network.\n",
    "transitions = np.zeros((n_items, n_items))\n",
    "for _ in range(n_samples):\n",
    "    src = np.random.choice(n_items)\n",
    "    neighbors = list(graph.successors(src))\n",
    "    if len(neighbors) == 0:\n",
    "        continue\n",
    "    dst = choix.compare(neighbors, params)\n",
    "    transitions[src, dst] += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The network looks as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd1hTZ/sH8G8gYcsSUAQUcACiIMOFMlTcYt211dbWqq21KtoBbhTrqNU6XlcddSstdb6uqgiOCoLIUBlSQWUpoGwCJDm/P/yZ15QhI+Fk3J/r4lKTk5M7GPLlGed5OAzDMCCEEEJUhBrbBRBCCCEtiYKPEEKISqHgI4QQolIo+AghhKgUCj5CCCEqhYKPEEKISqHgI4QQolIo+AghhKgUCj5CCCEqhYKPEEKISqHgI4QQolIo+AghhKgUCj5CCCEqhYKPEEKISqHgI4QQolIo+AghhKgUCj5CCCEqhYKPEEKISqHgI4QQolIo+AghhKgUCj5CCCEqhYKPEEKISuGyXYAs5JdWIvReJpJzi1HMF0Bfiwv7tvqY6GaJ1nqabJdHCCGERRyGYRi2i5CW+OeF2B6ehojUPABApUAkvk+LqwYGgI+dKb727gRnK0OWqiSEEMImpQm+I5EZ+PFCMvgCIep7RRwOoMVVx5IR9pjax7rF6iOEECIflKKr803oJaGiWvTeYxkGqKgW4scLSQBA4UcIISpG4Vt88c8LMXlPJCqqheLbiu+dQ1niNVTlZUDXwRsmoxbU+lhtnjpCZvWBkyV1exJCiKpQ+Fmd28PTwBcIJW7j6rWGgceH0HMaXO9j+QIhdoSnybI8Qgghckahgy+/tBIRqXk1xvR07Dyg06Uv1LT16308wwDXU/JQUFopwyoJIYTIE4UOvtB7mc0+BwdAaGzzz0MIIUQxKHTwJecWS1yy0BR8gQjJOSVSqogQQoi8U+jgK+YLpHSeaqmchxBCiPxT6ODT15LO1Rj6WjypnIcQQoj8U+jgs2+rD01uzZfAiIRgBFWASAgwIjCCKjAiYS1neLOii715K1mXSgghRE4o9HV8+aWV6Lc+rMY4X+HNoyi6fVziNoN+H8HQc0qNc2hy1fB3wEBaw5MQQlSEQgcfAMw6HIMrSS/qXaasLhwOMLRrG+ya6i79wgghhMglhe7qBIA5Pp2gxVVv0mM11Dn42qeTlCsihBAizxQ++JytDLFkhD20eY17KVyIwP/7GNSLsmRUGSGEEHmk8MEHvFloeskIB2jz1MHh1H8sh/Nmjc6gD7pjzfThGDhwIG7cuNEyhRJCCGGdwo/xvSshsxA7wtNwPSUPHLy5OP2tt/vxDbAzxdc+ncQLU1+9ehUff/wxduzYgQkTJrBTOCGEkBajVMH3VkFpJUJjM5GcU4JifjX0tXiwN2+FCa6178AeFxeHUaNG4YcffsC8efNYqJgQQkhLUcrga4qMjAwMHz4cfn5+WLduHdTUlKIXmBBCyL9Q8L3j1atX8PPzg7W1NX777TdoaGiwXRIhhBApo2bNO4yNjXH16lWUl5djxIgRKCoqYrskQgghUkbB9y/a2toIDQ2FnZ0dvLy8kJ2dzXZJhBBCpIiCrxbq6ur4z3/+g48++ggeHh5ISkpiuyRCCCFSIp3tDZQQh8NBYGAg2rVrBx8fH/z555/o378/22URQghpJprc0gB//fUXpkyZgt27d2PcuHFsl0MIIaQZqMXXAEOGDMHly5fh5+eH7OxsfPPNN2yXRAghpImoxdcI6enpGDZsGMaNG4c1a9aA87710QghhMgdCr5Gys/Ph5+fHzp37oy9e/fStX6EEKJgaFZnI5mYmODatWsoKirCyJEjUVxczHZJhBBCGoGCrwl0dHTw559/omPHjvD29kZOTg7bJRFCCGkgCr4m4nK52LlzJyZMmAAPDw8kJyezXRIhhJAGoFmdzcDhcLBkyRJYWFjAx8cHJ0+ehIeHB9tlEUIIqQdNbpGSixcv4tNPP8WePXswZswYtsshhBBSB2rxScnw4cNx8eJFjB49Gjk5OZg9ezbbJRFCCKkFtfik7MmTJxg2bBgmTpyI1atX07V+hBAiZyj4ZCAvLw+jRo2Cg4MD9uzZAx6Px3ZJhBBC/h/N6pQBU1NThIWFoaCgAH5+figpKWG7JEIIIf+Pgk9GdHV1cerUKbRv3x4+Pj7Izc1luyRCCCGg4JMpLpeL3bt344MPPoCHhwdSU1PZLokQQlQezeqUMQ6Hg+XLl8PCwgJeXl44ffo0+vTpw3ZZhBCismhySws6f/48PvvsM+zfvx9+fn5sl0MIISqJujpb0MiRI3H+/HnMmjULu3fvZrscQghRSdTiY0FaWhqGDRuGjz/+GCtXrqRr/QghpAVR8LHk5cuXGDlyJLp3747du3fTtX6EENJCKPhYVFpaikmTJgEAfv/9d+jp6bFcESGEKD8a42ORnp4ezpw5A3NzcwwYMAAvX75kuyRCCFF6FHws4/F42Lt3L4YPHw4PDw+kpaWxXRIhhCg1uo5PDnA4HKxatQqWlpbw9PTEmTNn0KtXL7bLIoQQpURjfHLm3LlzmD59Og4cOICRI0eyXQ4hhCgd6uqUM35+fjh37hy++OIL7N27l+1yCCFE6VCLT06lpqZi2LBhmDZtGpYvX07X+hFCiJRQ8Mmx3NxcjBw5Eq6urti5cye4XBqSJYSQ5qLgk3MlJSWYOHEiuFwuQkJCoKury3ZJhBCi0GiMT861atUK586dg4mJCQYOHIi8vDy2SyKEEIVGwacAeDwefvvtNwwePBgeHh74559/2C6JEEIUFg0aKQgOh4PVq1fDwsICnp6eOHv2LNzd3dkuixBCFA6N8Smg06dPY+bMmTh06BCGDx/OdjmEEKJQqKtTAY0ZMwZnzpzB559/jt9++43tcgghRKFQi0+BpaSkYNiwYfjiiy+wZMkSutaPEEIagIJPweXk5GDEiBHo3bs3/vOf/9C1foQQ8h4UfEqguLgY48ePh46ODo4fPw4dHR22SyKEELlFY3xKQF9fH+fPn4e+vj4GDRqE/Px8tksihBC5RcGnJDQ0NHDo0CH4+PigX79+SE9PZ7skQgiRSzQgpEQ4HA7Wrl0LCwsL9O/fH+fOnYOrqyvbZRFCiFyhMT4ldfLkSXz55Zc4evQohgwZwnY5hBAiN6irU0mNGzcOp06dwieffIJDhw6xXQ4hhMgNavEpuaSkJAwfPhyzZs3CokWL6Fo/JZNfWonQe5lIzi1GMV8AfS0u7NvqY6KbJVrrabJdHiFyiYJPBWRnZ2P48OHo378/tm7dCnV1dbZLIs0U/7wQ28PTEJH6ZreOSoFIfJ8WVw0MAB87U3zt3QnOVoYsVUmIfKLgUxFFRUUYN24c9PX1cezYMWhra7NdEmmiI5EZ+PFCMvgCIer76eVwAC2uOpaMsMfUPtYtVh8h8o7G+FSEgYEBLl68CB0dHfj6+qKgoIDtkkgTvAm9JFRU1x96AMAwQEW1ED9eSMKRyIwWqY8QRUAtPhUjEokQGBiIs2fP4tKlS7C2tma7JNJA8c8LMXlPJCqqhTXuK3sUgcLbxyEszoO6rhFaj/SHllU38f3aPHWEzOoDJ0vq9iSEruNTMWpqavjpp59gYWGBfv364fz58+jRowfbZZEG2B6eBr6gZuhVpN/H6/ADMP0gABrtukBY+qrGMXyBEDvC07BrKu3hSAh1daqo+fPnY8uWLRgyZAiuXr3KdjnkPfJLKxGRmldr92bRraMw6PcRNC3sweGogdvKBNxWJhLHMAxwPSUPBaWVLVQxIfKLgk+FTZgwAaGhoZgyZQqOHDnCdjmkHqH3Mmu9nREJUZmTBlF5EbJ2zUTm9ml49ddOiKprBhwHQGhs7echRJVQ8Kk4Ly8vhIWFYcmSJfjpp59AQ77yKTm3WOKShbeEZYWASIDylNtoM3U9zD/fiqoXT1D0d0iNY/kCEZJzSlqiXELkGgUfgaOjI27fvo3Dhw9j/vz5EAprjiORlldWVob79+/jxIkTeP6i9lm4HN6bi9RbufmBq2cMdR0DtOo5BhX/xNR6fDG/Wmb1EqIoaHILAQBYWlri5s2bGDt2LCZNmoQjR47QtX4tgGEYZGZmIiUlBcnJyUhJSRH/PS8vDzY2NjAwMICG90yA06bG49W19KD+r/G8+uhr8aRZPiEKiYKPiBkaGuLSpUuYNm0ahgwZgjNnzsDY2JjtspRCeXk5UlNTxeH29s/U1FS0atUK9vb2sLOzg62tLdq1a4eOHTsiOjoaKSkp6NChAzqYGyA3X63W7k697r4oufdfaNu6AepclMScgU6nnjWO0+Kqwd68VUu8XELkGl3HR2oQiUT4/vvvcfHiRVy6dAmampoIDAzE7t27oaGhwXZ5cothGGRlZdUIt5SUFLx8+RKdOnWCnZ2dOOTs7OzQuXNnpKen4+rVq7h69Sru3LmDbt26wdfXF76+vujTpw80NTWRX1qJfuvDag0+RijAq6u/ouxRBDhcHnTtPWE04HNwuJL/V5pcNfwdMJDW8CQqj4KP1OmXX37Bzz//DA0NDTx//hwhISEYP34822Wx7m3r7d1wS05OFrfe/h1u9vb26NChg3iN1CdPnoiDLiwsDKampuKg8/b2hqFh7ReZzzocgytJL967YkttGJEIFY8jUf7XFhgbG0NPTw9t27bFyZMn63w+QpQVBR+pU3V1Nbp3746UlBQAQJ8+fXDnzh2JY5R1d4C3rbd/h9vb1lvHjh3F4fZuyBkYGNQ4V35+PsLCwsRhV15eLg66QYMGwcrKqkE11bdyy/tocdVQ8Psy5KfeE99maWmJJ0+egMejcT+iWij4SJ2WLVuGNWvWQCR6073G4/GQnJwMW1tbpdkd4N3W27+7J/X09GoNN2tr63p3uKioqMCtW7fEQff48WN4eXmJw87R0bHJ20P9b63Oml2eddHmqWHJCAf0bws4ODiAz+cDAFxcXHDx4kW0aVNz0gwhyoyCj9QpKysLISEhOHnyJKKjo1FVVQU/Pz9MWrxVoXYHqK319vbvb1tv73ZL2tvbo0uXLg3uAhQKhYiNjRUHXVRUFHr06CEOul69ekl1bLQ5uzPs2LED8+bNg62tLfLz81FaWop58+Zh9erV0NLSklqNhMgzCj7SIJWVlTh37hzu5HNxIVuzSS0OWYdfXa231NRU6Orq1hh3a0jrrTYMwyAtLU0cdNevX0e7du3EQefl5QV9fX0Zvco3EjILsSM8DddT8sDBm4vT33rb4h5gZ4qvfTpJLEzNMAzmzJmDOXPmoHPnzli2bBk2b94MbW1tbN26FZ988gltVkyUHgUfabD6xphyjwaiMjsFHLU3IaLeqjUsZu0W31/X7gBFRUXQ19dv8IdtQ1tv73ZRNqb1VpeXL1/i2rVr4rATCAQS43Tt2rVr1vmbqqC0EqGxmUjOKUExvxr6WjzYm7fCBNeGj7FmZ2dj2rRpiIiIgI2NDQ4dOoTevXvLuHJC2EPBRxqsvlmFuUcDodttAFo5D631sRwOMLRrG/HuAAzDYO/evZg7dy7OnDmDoUMlH9fQsbfmtt7qUlZWhps3b4qDLiMjA97e3uKws7e3V7qW0Y0bNzB16lS8ePECvr6+2LVrV4Mn3hCiSOgCdtIg9e0O0BASuwNUlmLq1Km4ceMGBAIBTp48ibS0tFpnTr4NtmHDhsHf318qrbfaCAQCxMTEiIMuJiYGbm5u8PX1xc6dO9GzZ09wucr94+Ll5YX09HRs27YNS5cuRefOnTF37lysWLECenp6bJdHiNRQi480yK6If/DL1dRaL6AG3rT4qvOfAQB4xhYw9PoEWh2cJI7R4qphlDUH277+AFVVVeLbzczMMGbMGJm13mrDMAxSUlLEQRceHo727dtj8ODB8PX1haenp0p/2BcUFMDf3x+hoaHQ0tLCxo0b8dlnn0FNjZb3JYqPgo80iH/IfZyOy67z/srsFPBaW4GjzkNZ0g28urIL5p9vBc/IXOK4gTZ6yDixEjExMRCJRKioqICjoyMSEhJk/RKQm5srMU4HQBx0AwcORNu2bWVeg6KJi4vDZ599hrS0NFhYWODXX3+Ft7c322UR0izK3XdDpKaYL6j3fs12duK/63UfhLJHEaj4JwY8dz+J4ziaOggLC4NIJEJMTAx2796NgoLadx5orpKSEty4cUMcdJmZmRgwYAB8fX2xaNEidO7cWenG6aStR48euH//Po4ePYr58+dj5MiR8PLywrZt29CxY0e2yyOkSSj4SIPoazXyrcLhAKjZmaDDBQ4dOoTw8HBcvHgReXl5cHBwkEqN1dXVuHv3rjjo7t+/j169esHX1xf79u2Dq6ur0o/TyQKHw8HUqVPxwQcfYOXKldi5cyecnJwwa9YsrFixgpY8IwqHujpJg9Q3xifil6IyOwVa7bsDaupvujov/Qfmn20Br7Wl+Dgtrhq64xlCV38l8XhfX18cP34cJiYN314HeDNO9+jRI3HQ3bhxA7a2tuKZl56entDR0WnaCyZ1Sk1NxezZs3Hv3pvlz9asWYNZs2bRLxVEYVDwkQapb3cAYXkRXv4ehOpXmQBHDbzWljD0nAptGxeJ497uDnDlv6cwffp0VFRUgMvlomPHjsjNzYWhoSHc3Nzg7u4ONzc3uLm5oXXr1hLnyMrKkhin09DQkBinMzU1len3gbzBMAz++9//4uuvvwafz4eBgQG2bduG4cOHs10aIe9FwUcarDm7A/z7Or5Hjx5h2LBhyMrKQkxMDJydnfHPP//g3r17iImJwb179xAbGwtDQ0NYWlqKL1wvKSnBoEGDxK06W1tbGqdjEZ/Px88//4yffvoJPB4Prq6u2Lx5MxwdHdkujZA6qQcFBQWxXQRRDO2NdXA6LhsCUeOTT5unjp/GO6GN/pv1IE1NTTFjxgwYGhrigw8+gLq6Olq3bo0uXbpAV1cXQqEQpaWlSE9Ph46ODnR0dMDj8ZCXl4fi4mLw+XwUFRWhsrISRkZGtFs8S7hcLry8vDB16lSkpqbi9u3b2Lt3L549e4bevXtDV1eX7RIJqYFafKRRmrM7QG1rdTIMg8TERHHX5a1bt9ClSxdxi65fv34SoSYSifD48WNxqzAmJgZxcXEwMzOT6CZ1dXWlSRcsuH79OubMmYPS0lIUFxdjyZIlmDdvHjQ1FXeLKqJ8KPhIozVndwAAePbsmcQ4XatWrcRBN2DAgBrjeu8jFAqRmpoq0U0aFxeHtm3bioPQ3d0drq6uMl88mrxZBWfHjh0ICgqCkZERRCIRfv75Z4wbN466pYlcoOAjTfK+3QGqBQK0EeVj19yxsNJlEB4eLg66V69eicfpBg0aBBsbG6nXJxQKkZKSItEyjI+Ph4WFhUTL0MXFhcJQRl6+fInFixfj1KlT0NHRgY2NDX755Re4ubmxXRpRcRR8pFlq2x2gk6k2NLPjMHfmZ3B0dERKSgo8PDzEsy+dnJxYWfpKIBAgOTlZomWYkJAAS0tLiZahi4uLSi9XJm3R0dGYM2cOCgoKUFRUhFGjRmHNmjWs7WhBCAUfaTaRSISEhARcuXIFV69exd9//42uXbuipKQEAwcOxM8//yy3m5wKBAIkJSXVCMP27dtLhGGPHj0oDJtBJBLh4MGDWLRoEdq2bYtnz57B398f3333HV1rSVocBR9pkoyMDHHX5bVr12BsbCwep/Px8YGRkREiIyPx8ccf4/HjxzJdcFraqqurkZSUJNFN+uDBA1hbW0t0k/bo0YNmLTZSUVERgoKCcPjwYbRv3x4vX77E2rVrMWXKFFoAm7QYCj7SIAUFBbh+/bo47EpKSsRjdL6+vmjfvn2tj+vTpw8CAwMxZsyYFq5Yuqqrq/Hw4UOJluHDhw9hY2Mj0TJ0dnamFkwDPHr0CPPmzUN6ejo0NDTQqlUrbNq0Cf3792e7NKICKPhIrSoqKnD79m1x0KWmpqJ///7iVl337t0bNEMvJCQEO3bsQERERAtU3bKqqqrw8OFDiZbho0eP0LFjxxphSNcZ1sQwDE6ePIlvv/1W3P3Zv39/rF+/XiYTngh5i4KPAHgzCzIuLk48ThcZGQknJydx0PXp0wcaGhqNPm91dTU6duyI06dPw9XVVQaVy5fKyko8ePBAomWYlJSEzp07S3STOjs7y+24Z0srLy/H+vXrsX37dvFuEDNnzsTixYtpxi2RCQo+FcUwDJ48eSJu0YWFhaFNmzbioPP29oaBgYFUnuunn37CgwcPcOjQIamcT9FUVlYiMTFRIgyTk5PRpUsXiZZh9+7dVToM09PTsXDhQsTFxaFjx4548OABVq5ciS+++IIWwCZSRcGnQvLy8hAWFiYOu8rKSnHQDRo0CBYWFjJ53tevX8PW1haPHj2Cubn5+x+gAvh8PhITEyW6SVNTU2Fvby9eoPttGKraqid//fUX5s2bBxMTE1RVVaGiogKbNm3C4MGD2S6NKAkKPiVWXl6OmzdvioPuyZMn8PLygq+vLwYPHgwHB4cWW0ljzpw5MDY2RnBwcIs8nyKqqKhAQkKCRMvw8ePHcHBwkOgm7d69e5O6nRVJVVUVtm7dinXr1sHLywvx8fGwt7fHxo0bYW9vz3Z5RMFR8CkRoVCIe/fuicfpoqOj4eLiIm7V9erVCzwej5XaUlJS4OnpiadPn9JEj0aoqKhAfHy8RMvwn3/+QdeuXSW6SR0dHZUyDHNychAQEICwsDD4+Pjg0qVL+OijjxAUFNTope0IeYuCT4ExDIPHjx+LW3TXr1+HpaWlOOi8vLzQqlUrtssUGzVqFMaMGYMZM2awXYpCKy8vrxGG6enpcHR0lGgZOjo6svaLjrTdvn0bc+fOhaamJjp06IBr165h8eLFmDNnjlIGPpEtCj4F8+LFC4kFnkUikcRGrPI8hnb16lX4+/sjMTGRFiuWsrKyMsTFxUl0k2ZkZKBbt24SYdi1a1eFDUOhUIi9e/di+fLlGDhwIPLz8/H06VNs2LABo0ePpvcUaTAKPjlXWlqKGzduiIPu2bNn8PHxEYddly5dFOYHnmEYODk50USFFlJaWoq4uDiJluGzZ8/QvXt3iW5SBwcHhZo1+erVKyxbtgyhoaH48MMPce3aNbRp0wabNm1Cjx492C6PKAAKPjkjEAgQHR2Nq1ev4sqVK4iNjYW7u7u4+9Ld3V2hPqT+bd++fTh58iTOnz/PdikqqaSkBPfv35doGWZmZsLJyUmiZWhvby/377P4+HjMnTtXvIrQ4cOHMXLkSKxevVquez4I+yj4WMYwDJKTk8UtuoiICFhbW4uDztPTU6nWg6yoqIC1tTUiIiJodp6cKC4uxv379yVahtnZ2XB2dq4RhvK25irDMDh+/Dh++OEH9OvXD61bt8bvv/+OBQsWYOHChTSRitSKgo8F2dnZEuN0XC5XYpzOzMyM7RJlavny5cjPz8eOHTvYLoXUoaioCLGxsbh37544DHNzc+Hs7CzRTdqlSxe5CMPS0lKsXr0ae/fuxYwZM/D48WPExMRg3bp1mDx5ssIMB5CWQcHXAoqLixERESEOupycHAwYMEB8PV3Hjh1V6gczJycHXbt2xT///ANjY2O2yyENVFhYKA7Dt63DFy9ewMXFRaJl2KVLF9Z2WkhNTYW/vz+ePHmCmTNn4vjx4+DxeNi0aRP69u3LSk1E/lDwyUB1dTWioqLEQRcXF4fevXuLuy9dXV3l4rdkNk2bNg1du3ZFQEAA26WQZnj9+jViY2Mluknz8/Ph4uIiDkI3Nzd07ty5xcKQYRicP38e/v7+cHR0hKenJzZv3gxPT0+sW7cOHTp0aJE6iPyi4JMChmHw8OFDcdDdvHkTnTp1Egddv379aKuaf7l//z5Gjx6NJ0+eKOz0elK7V69eibtI34bhq1evJMLQ3d0dHTt2lGkY8vl8bNq0CZs2bcKMGTOgrq6O3bt348svvxQviabqv4CqKgq+JsrMzBQH3dWrV6GrqysOugEDBsDExITtEuWej48PvvrqK0yePJntUoiMFRQUSAThvXv3UFhYCFdXV4luUll0+z9//hzff/897ty5g8WLF+PWrVs4ceIE+vTpg/Dw8DrDL7+0EqH3MpGcW4xivgD6WlzYt9XHRDdLtNZTrfVTlY3CBB/bb8KioiKJjVjz8/MxcOBADB48GIMGDYKtra3Ma1A2p0+fxtq1axEZGalSY5zkjby8vBrdpCUlJXB1dZVoGdrY2Ejl/REeHo65c+eiuroaGRkZqK6uho2NDX799VcMHDhQfFz880JsD09DRGoeAKBSIBLfp8VVAwPAx84UX3t3grOVYbPrIi1P7oOPrTdhZWUlIiMjxUH34MED9O3bV9yq69GjB2sD+MpCKBSiS5cuOHLkCE08IACAly9f1mgZlpWViccK34ahtbV1k8IwPz8fVlZW4PP5AAAjIyMYGBjAyckJGzZswN1XGvjxQjL4AiHq+2TkcAAtrjqWjLDH1D7WTXy1hC1yHXxHIjNa7E0oEomQmJgoDrrbt2/D3t5eHHQeHh4qvVearGzZsgW3b9/G77//znYpRE69ePFCIghjYmLA5/MlgtDNzQ0dOnR4bxguXboU69evB4/HA5/PB8MwGDVqFLy8vLDpbDR0+k2FkNPwcT9tnhqWjHCg8FMwcht8b0IvCRXVovcf/P9qexMWFxfjwIEDmDt3bo0fiqdPn4qD7tq1azA0NBQHnY+PD021bwElJSWwtrZGbGwszbYjDZaTk1OjZVhVVSURhO7u7rCyspL4uX/27Bnu3r2L169f49WrV4iPj4empib8V/2CD3+9A76g5udN9assZO/7Brr2/WDi912N+7V56giZ1QdOltTtqSjkMvjinxdi8p5IVFQLxbcxgmoU/LUD/Iw4iPil4Bqaw8j7U2h3dJd47LtvwqysLHh7e+PJkydISkqCqampxDhdUVGRxEas9MHLjoULF0JdXR0bNmxguxSiwLKzsyXCMCYmBkKhUCII3dzcYGlpWeOX4FmHY3Al6UWtPUsvTiwDI6gE18Cs1uDjcIChXdtg11T3mg8mckkug6+2N6Goio/iqD+h190X6gamqPgnBvlnN6Dd9P+Aa9hGfNzbN+FcF214e3NHeNsAACAASURBVHvj9evXUFNTg7m5OV6/fo3+/fuLw6579+40TicH0tPT4e7ujqdPn0JPT4/tcoiSYBhGHIbvdpMCkOgmdXTtg5G/xkrMH3ir7FEEylPvgNfaCoLCnFqDDwA0uWr4O2AgzfZUEHIXfPmllei3PqzWN+G/Ze/7Bgb9PoKufT+J23lqQMaWqRCUFYpvc3Z2RlRUFDQ16Y0pj8aPH48BAwbgm2++YbsUosQYhkFWVpZEEGq7+uGhmk2NzxxRZTlyDvijzUc/ojT+r3qDT4urhgWDu+BLr44t8TJIM8nd8uuh9zIbdJyw7DWqX2VBw7R9jfs4HA76ffIdmEd/4f79+ygpKRH35RP5tGDBAnz++ef4+uuvqRVeD4ZhxF8ikUji39K+TdbnZ7u29u3bw8rKCpEcO1Q+Ka3xvS68cRh6zkPA1Td97/8LXyBCck6JLP7LiQzIXfAl5xa/t7XHCAXIP/sz9LoPAq+1VY37q4QMTDs74/GdPyAQCAAA6urqePDggdz9ALP9wy8vt4lEIhQUFMDLywvm5uYq831qzOPepaamBg6HI/ElzdtkfX55qq3AvpbPkBdPwH8aD/PPtzT4s6uYX93gYwm75C74ivmCeu9nGBHy/7sRUOfCePBXdR53MzIGL+Ljxf8WCoX48MMP5eoHWN4+ENTV1Vl9DaampggPD8f8+fPl5nsij/9fRLr8Q+4jPj9b4jb+s0QIil4gc8fnAACmig8wIuTkz68zDPW1aOk9RSF3waevVXdJDMOg4MJWCMsKYTYxCBz1uo+dNMYPth6GWLVqFSorK1FVVYWHDx/KomQiJePHj4e1tTW6dOkCZ2dntsshKsK+rT40ubkSPU16PYZC18FL/O/iuychKHoB46Fzaj0HRyRAW21hrfcR+SN3gylv3oS1l/Xq8nZUFzyH2YTlUOPVPV6nxVWDo6Uh5s+fj4SEBHzzzTfo3bu3rEomUqKhoYE5c+Zgy5aGdy8R0lwT3Cxr3KbG04K6npH4i8PTAoerAXUdg1rPwVFTw4avJ2D58uUoLa05Xkjki8LM6hQUvUTWzumAOg8ctf+trGA8bA70HAdIHMsIqpC7ewaE5UXgcrlgGAZ9+/ZFeHh4S7wE0gz5+fno3LkzkpOT0aZNm/c/gBApqO86vvcSicDLS8JKX0tcuHAB169fx+rVqzFt2jSaqCWn5O5/xURPE95dTPHvoQyugRk6BP4XHb4/hfbfhoq//h16HA7QzZgDYXkRhEIhKisrIRAI0KtXrxZ8FaSpTExMMHHiROzatYvtUogKmePTCVrcpm1RpK3JxVdetlixYgVevHiBNWvWYO/evXB3d0dERISUKyXSIHfBBzTvTajFVcfaT3xw9+5d6OrqAngzNrhhwwbo6+tjypQpSE9Pl2a5RMr8/f2xc+dO8ULChMias5UhloywhzavcR+Jb5dJ/PazCXj06BHGjx+PRYsWwcrKCtOmTcO0adMwbtw4pKWlyahy0hRyGXzNexPaw8nSEK6urrh8+TJ4PB6GDx+O0tJSzJ8/Hzdu3ICtrS0sLCwQEBCA4uJiGb0K0lRdu3aFs7Mzjh8/znYpRIVM7WONJSMcoM1Tr9Hj9G8czpvlEd9dG5jH4+Grr77C48eP0b17dwQHB2PYsGFwcHBAnz598N1336GwsLD+E5MWIXdjfO+Sxu4M0dHRMDQ0ROfOncW3PX36FCtWrMDp06dRXFwMBwcHzJ07F7NmzaI+eTlx6dIlBAQEIC4ujqbwkxaVkFmIHeFpuPIwB2ocDqqZ/73/3m6FNsDOFF/7dKp3Yer8/HysWbMGBw8exLRp0/Dq1StcvHgRy5cvx5dffgkuV+4m1asMuQ4+4H9vwuspeeAAEqunN+ZNWJfbt29j1apVCA8Ph0gkgoeHB5YuXYrBgwdL70WQRhOJRHB0dMSOHTswYMCA9z+AEClz8/DGkC+Xga9lgmJ+NfS1eLA3b4UJro3b/Prp06dYvnw5Ll26hGnTpiEmJga5ubnYuHEjhg8fLsNXQOoi98H3VkFpJUJjM5GcU9KsN2FdRCIRjh07ho0bNyIhIQHa2toYMWIEgoODYWdnJ4VXQBpr9+7dOH/+PM6ePct2KUTF8Pl8GBsbIz8/Hzo6OlI5Z2JiIhYtWoQHDx5g3LhxOH/+PGxsbPDzzz+jW7duUnkO0jAKE3wtic/n46effsL+/fvx9OlTtGnTBlOmTMGyZctgaEh7brWU8vJydOjQAX///bdEVzUhshYZGYnZs2fj/v37Uj/3zZs3ERAQgJKSEnh4eODkyZOYMGECVq1aBVPT968LSpqPBrRqoaWlheXLlyMjIwOZmZnw8/PDb7/9BmNjY9jb22Pbtm3iNUCJ7Ojo6GDmzJnYunUr26UQFXP37l2ZXQLl6emJ27dvIzg4GDdv3kSXLl3w+vVrdO3aFRs2bEBlZaVMnpf8DwXfe1hYWGDPnj149eoVoqKiYGtri++//x5aWlro168fLly4wHaJSm3OnDk4cuQIzYYjLSoqKkqmqz1xOByMGTMGCQkJmD59Om7fvg0XFxdcvHgRDg4OCA0NBXXGyQ4FXyP07NkTFy5cQEVFBY4dO4bKykr4+flBV1cX48aNo7VAZcDCwgIjRozA3r172S6FqBBZtvjexeVy8cUXXyA1NRWDBw9GYmIi7O3tsXz5cnh7e4s3ziXSRcHXBBwOB5MmTUJMTAwqKiqwbNkyxMfHo1u3bjAzM8O8efOQn5/PdplKw9/fn7qXSYspKCjAixcv4ODg0GLPqa2tje+//x6pqalwcnLCixcvoKOjg1GjRmHatGnIyspqsVpUAQVfM2loaCAwMBD//PMPXrx4gfHjx+PYsWMwMzND586dsXHjRvrAbqaePXvCysoKp0+fZrsUogKio6Ph5uYGdfWmrR7VHEZGRli3bh3i4+NhZWUFgUCA58+fw8nJCStXrkRZWVmL16SMKPikyMzMDDt37kR+fj5iY2PRtWtXLFu2DFpaWujduzdOnTrFdokKy9/fH7/88gvbZRAVIOvxvYawtLTEnj17cPPmTRgaGkJTUxPnz5+HnZ0dDh06VGNjYtI4FHwy0qNHD5w5cwbl5eX4448/wOFwMHHiROjo6GD06NGIi4tju0SFMmbMGGRlZeHu3btsl0KUXEuN7zWEg4MDTp48iZMnT0JbWxs8Hg+rV69Gr169cPPmTbbLU1h0HV8LEggE2LJlC3bv3o20tDQYGxtj4sSJCAoKoi14GmDjxo24d+8ejh07xnYpREkxDANTU1PEx8fDwsKC7XIkMAyDixcvIjAwEHw+H0VFRfD09MRPP/0EW1tbtstTKBR8LMnPz8eqVatw4sQJ5OXlwcbGBjNnzsS3334LDQ0NtsuTS0VFRbCxsUFCQgIsLWtuHkpIcz158gSenp5yPZlEKBTi2LFjWLp0KbS1tZGbm4tZs2ZhyZIlMDCofaNcIom6OlliYmKCrVu34uXLl3jw4AFcXFywevVqaGtrw93dHSEhIdSP/y8GBgaYOnUqtm/fznYpREnJw/je+6irq+OTTz5BamoqvvrqK/B4PJw6dQqdOnXCrl27aDJdA1DwyQFHR0f8+eefKCsrw7lz56ClpYUpU6ZAR0cHI0aMQHR0NNslyo158+Zh7969KC8vZ7sUooTkaXzvfTQ1NeHv74+0tDRMnjwZ1dXV+PHHH9GtWzdcvnyZ7fLkGgWfnBkxYgRu3boFPp+PDRs24MmTJ+jduzeMjY0xY8YMZGdns10iqzp16gQPDw8cOnSI7VKIElKEFt+/GRgYIDg4GMnJyRg1ahSys7Px8ccfY8iQIUhKSmK7PLlEY3wKoLCwEMHBwTh27Bhyc3PRvn17TJ8+HQEBAdDS0mK7vBYXHh6O2bNn4+HDh7R/IpGa6upqGBoaIicnB/r6+myX02RpaWlYvHgxLl++DKFQiE8++QTBwcEwMTFhuzS5QZ8aCsDQ0BAbN25ETk4OkpOT0bdvX/z888/Q1dWFi4uLyl3X4+3tDU1NTerOIVKVkJAAGxsbhQ494E2vyO+//46wsDC4u7vj+PHjsLW1xcaNG1FVVcV2eXKBgk/B2NnZ4cSJEygpKcHly5ehr6+PL774Atra2hg8eDBu377Ndokyx+FwsGDBAmzevJntUogSuXv3rsJ1c9bHzc0N4eHh+OOPP2BlZYXg4GBYW1vj1KlTKr8ANgWfAvP19UVERAQqKyuxbds2ZGdnw9PTE4aGhpg2bRqePn3KdokyM3nyZCQkJNDC4ERqoqKiFGZiS2O8Xfx69+7d4HA4+PTTT9GzZ0+Z7DWoKCj4lICamhpmzZqFhw8foqioCLNnz8bVq1dhbW0NKysrLFmyBKWlpWyXKVWampqYPXs2tmzZwnYpREko4sSWhlJTU8OHH36I9PR0rFmzBo8fP0bfvn0xfvx4lZwwR5NblFh6ejqWLVuGc+fOoaSkBN26dcP8+fPx+eefK8WkkJcvX8LOzg6PHz+mgXvSLEVFRbCwsEBhYSG4XC7b5chcaWkp1q5di19++QUikQjz58/HihUroKOjA+DNKlPK/H1Q/E8/UicbGxscOXIERUVFuH79OkxNTTF79mxoaWlh4MCBCA8PZ7vEZjEzM8O4ceOwe/dutkshCi46OhouLi5K/WH/Lj09Pfz44494+vQpPvroI2zevBnm5ubYs2cPnjx5gtatW793Xdz80krsivgH/iH3Mf1gNPxD7mNXxD8oKJX/HeSpxadiRCIRDh48iM2bNyMxMRF6enoYNWoUgoOD0bFjR7bLa7TExEQMHToUGRkZtNQbabI1a9agoKAAGzduZLsUVmRkZOCrr77CtWvXwOPxwOfzYW5ujpSUFOjp6UkcG/+8ENvD0xCRmgcAqBT8b0a5FlcNDAAfO1N87d0JzlaGLfkyGoxafCpGTU0Nn3/+OeLj41FaWgp/f3/cunULnTp1Qrt27RAQEIDi4mK2y2yw7t27o2vXrvj999/ZLoUoMGUe32sIa2trXLp0CUePHgWfzwfDMMjJycEnn3wicdyRyAxM3hOJK0kvUCkQSYQeAPD//7a/Hr3A5D2ROBKZ0YKvouEo+FSYjo4OVq1ahWfPnuHp06cYOnQofv31VxgaGqJr167YuXOnQlwf+HavPuq8IE3BMIzSzuhsrL1794LD4QB48305ffo0fH19UVxcjCORGfjxQhIqqoV4348awwAV1UL8eCFJLsOPujpJDXfu3EFQUBDCw8MhEonQt29fLFu2DIMHD2a7tFqJRCLY29tj37598PT0ZLscomCePXuGnj17Ijc3V/yhr6rOnDmD5ORkVFRUoLy8HDExMSgpKUG/0R/jcrUDKqqFEsfnn/sZ/Ix4iKr5UNc1gn6f8WjlPFTiGG2eOkJm9YGTpfx0e1LwkTqJRCIcP34cGzduRHx8PLS1tTF8+HCsXr0adnZ2bJcnYfv27bh27RpOnjzJdilEwYSGhuLgwYM4d+4c26XIrVmHY3Al6UWNll5V3lPwjNqBw+WhuuA5co8tgtnEIGi27SQ+hsMBhnZtg11T3Vu46rpRVyepk5qaGqZMmYLY2FiUlZUhMDAQ0dHRsLe3R9u2bbFw4UIUFhayXSYAYNq0abhx4waePHnCdilEwaj6+N775JdWIiI1r9buTQ3TDuBwef//Lw444EDwOkfiGIYBrqfkydVsTwo+0iBaWlpYunQpMjIykJWVBT8/Pxw8eBDGxsaws7PDli1bWN0HTE9PD9OnT8e2bdtYq4EoJkXaikiWAgMDMWrUKNy8eVNivDz0Xma9jyu4vAPPfh6P7D1fQV3PGNoda7bsOABCY+s/T0uirk7SLDExMVixYgXCwsJQXV2NXr16YfHixRg1alSL1/Ls2TP06NEDGRkZaNWqlcqP15D3EwgEMDIywrNnz2BkZMR2OayaPHkyQkJCoKurCxMTE8ycORPz58/H0vOPcTqu/tVdGJEQlVnJ4D9LhEGfCeCo17wecmwPC/zyYQ9Zld8oFHxEakJDQ7F+/XrExsZCS0sLQ4YMQXBwMLp169ZiNYwfPx7Am+6rmTNnYsWKFS323ETxJCQkYOLEiUhJSWG7FKkSCAR4/fo1CgsL8fr16/d+FRYWIi0tDSUlJRLnWbNmDR6380VY8ssGPW/Bpf+AZ9Ie+u6ja9w3yN4M+6b1lMrray7VWKaAtIgJEyZgwoQJqKqqwqZNm7Bnzx44OTnBxMQEH374IVasWCGzpcVEIhECAwNx4cIFVFZWgmEYCIXC9z+QqDR5Ht+rqqpqUGDVdntFRQUMDAxgaGgIIyOjGl8mJibo3Lmz+N+GhoY4duwYNm/eDB0dHXTq1An79++Hq6sr/EMasZi1SFRjjO8tfS1erbezgYKPSJ2GhgYCAwMRGBiIly9fIigoCCdOnMD27dtha2uLL7/8EvPnz5fqSisikQiXL18GwzBgGAbq6uoKv68akT1Zj+9VVFQ0KrDe/Xq7MW5twWVkZIR27dqha9eutd7XqlWrRq/H++DBAxgZGWHr1q34+OOPxUMF9m31ocnNrXGxurCsEPyn8dDu1Ascrgb4GXEoS4qAid/3Nc6txVWDvXmrpn8jpYy6OkmLSUhIwPLly/HXX3+hsrIS7u7u+OGHH8Tdk83F5/MxdepUnDt3DtXV1di5cye+/PJLqZybKCcnJyfs27cPPXvW3gXHMAzKysoaHFb/PoZhmFqDqb5Ae/ulq6vbouPUIpEIIpGoxnql+aWV6Lc+rGbwlRch79RaVL1MBxgRuAZmaOXmh1Y9htU4tyZXDX8HDERrPU2ZvoaGouAjrDhz5gzWrVuH6Oho8Hg8+Pr6YtWqVXBxcWnWeRmGwcqVK7Fy5UqsXbsWgYGBUqqYKCqGYVBSUlIjpHJycrBgwQIsXLiw1vvfBhmXy21wWP37GG1tbaWYZPX5/jsIf/wKTQkLebyOj4KPsEogEGDr1q3YtWsX0tLSYGRkhIkTJyIoKAht27Zt8nkPHDgAb29vtDJth9B7mUjOLUYxXwB9LS7s2+pjopul3Pz2Sd5PJBKhqKioSV2GRUVF0NLSqhFSlZWViIuLw9dff11vqGlqqvb7JDExEX6ffQPe0O9RzTQ+xGnlFkLq8erVK6xcuRInTpzAy5cvYW1tjRkzZuDbb7+FlpaWxLEFBQVo3bp1nedShhXklY1AIBAHVWO7DYuLi6Gnp9ekLkNDQ0PweDUnVmzYsAGZmZm0mXE9zp49iy+++AJbtmyByNbj/9fqbPj6vdo8NSwZ4YCpfaxlV2QTUPARufTo0SMsX74cly5dQkVFBVxcXPDdd99h0qRJyMrKgo2NDXbu3ImZM2fWeOybxXSTwRfUv5guhwNocdWxZIS93P1gyqvq6uomzzQsKyuDgYFBk7oNDQwMpL5X3oQJEzB27FhMmTJFqudVBgzDYMOGDdiyZQtOnjwpnvmqLD9bFHxE7l24cAFr1qxBZGQkuFwu2rdvj/T0dPB4POzfvx+TJ08WH/u/FeQV/7dSWeHz+U3qMiwsLERlZSUMDQ3fG1y13a+vr9/omYay1L59e4SFhaFTp07vP1iFVFZWYtasWUhMTMTZs2dhaWkpcX9CZiF2hKfhekoeOHizFdFbb3tTBtiZ4mufTnLVvfkuCj6iMAQCAXbu3ImFCxeKl0dTV1fHnj173uwx+LwQk/dESqwg/2zjBIlzMIIqtHIZAeMhX0ncLo/jEHVhGAbl5eVNmmX4+vVrCIXCRofW2y89PT2lmKyRk5ODbt26IT8/Xylej7S8ePECY8eORbt27XDw4EHo6urWeWxBaSVCYzORnFOCYn419LV4sDdvhQmu8j9+TtfxEYXB5XLRv39/CIVCaGlpQSQSoaqqCtOnT8eVK1egM9QffIHkRevtvw0V/11UxUfmtqnQse9f49x8gRA7wtNabOYZwzAoLS1t0jVehYWFUFNTqzewbG1t6ww1HR0dlf+wf7v/nqp/H94VHx+PDz74AJ9++imCgoLe2zpvraeJL706tlB10kXBRxSKsbExAgIC0KFDB1hZWaF9+/bQ09ND+J0YrEsqqHfcoTzlNtR1DKBp5VjjvndXkG/ob6sikQjFxcUN7ip897jCwkJoamrW28qys7Or8/5/T/YhjUMLU0s6c+YMZsyYgW3btkkMHSgrCj6iUDp06IC1a9fWuP3yMxGQlFrvY0sTr0G328A6f8tnGAZrT4Shl35Jg1pfxcXF0NXVrbdr0NLSstb7DQ0NpbpyDWmcqKgofPvtt2yXwTqGYbBu3Tps374d58+fV5lfBij4iFJIzi2usbLEuwRFL1H5/AFaj5hX5zFVQgaXoh7iaeFdiYCytrauNdRkMdOQyJ5IJEJMTEydq7WoCj6fj5kzZyIpKQlRUVGwsLBgu6QWQz+1RCkU8+vfC7D0QRg0LbuCZ1j/RfG9+/tg37Saaw0S5ZGcnAwTExOYmpqyXQprcnNzMXbsWFhZWeHGjRvQ0dFhu6QWJT9ziwlpBn2t+n+HK3sQBr1uAxtwHvlZQZ7IhqqP78XFxaF3794YNmwYQkJCVC70AGrxESVR1wryAMDPTIKwtKDW2ZzvkrcV5IlsyPNWRLJ26tQpzJo1C9u3b8ekSZPYLoc11OIjSmGCm2Wd95U9uAadLh5Q06z/N1sGwATXus9DlIMqtvgYhsGaNWswb948XLx4UaVDD6AL2IkSmXU4BleSXtR7SUNdOACq0qNhknQK5ubmMDIyQps2bbBy5UoYGBhIvVbCjoqKCrRu3RoFBQXQ1tZmu5wWUVFRgRkzZiA1NRVnzpxBu3bt2C6JddTVSZTGHJ9OuPk4X2LllobS5KmhIuEiEpISkJCQAADQ0dFBUFCQlKskbLp//z4cHBxUJvRyc3MxZswYWFtbIyIiQiXH82pDXZ1EaThbGWLJCHto8xr3ttbmqWHpCAdcCdkrcWH44MGD6YNCybxdsUUV3L9/H7169cLIkSNx/Phxei+/g4KPKJWpfayxZIQDtHnqeN9qVBzOmzU63y5Q3b17dwQEBIDH46FNmzaoqqpCt27dcO7cOdCIgHK4e/euSkxs+fPPPzFkyBBs2rQJy5Yto6XZ/oXG+IhSauoK8gKBAL6+vli8eDGGDBmCy5cvY+HChTA3N8emTZvg5OTU8i+GSI2trS3Onz8PBwcHtkuRCYZh8OOPP2L37t04ffo03Nzc2C5JLlHwEaUmjRXkBQIBfv31V6xcuRJjx47FqlWrYGZmJuPKibTl5eWhc+fOePXqlVxtjyQtFRUV+OKLL5CWloYzZ87A3Nyc7ZLkFgUfIQ30+vVrBAcH4/DhwwgICMDcuXOhqSnf26+Q/zl//jx++eUXXL16le1SpC47OxtjxoxBx44dsX//fpWZvNNUyvdrDyEyYmRkhE2bNuH27du4efMmHB0dcerUKRr/UxDKOr4XGxuL3r17Y/To0Th27BiFXgNQ8BHSSF26dMGZM2ewa9cuLF++HAMHDkRcXBzbZZH3UMYZnaGhoRg6dCg2b96MpUuX0iSWBqKuTkKaQSAQYN++fVixYgX8/PywevVqtGnThu2yyL8wDIPWrVvj4cOHSjH2xTAMgoODsW/fPpw+fRouLi5sl6RQqMVHSDNwuVx8+eWXSElJgaGhIRwdHbF+/Xrw+Xy2SyPvSEtLg56enlKEXkVFBT766CNcuHABUVFRFHpNQMFHiBQYGBhgw4YNiIyMRGRkJLp27YrQ0FAa/5MTyrIwdXZ2Nry8vMDlchEeHo62bevfZovUjoKPECnq1KkTTp06hX379mH16tXw9vZGbGws22WpPGVYmDomJga9e/fGuHHjcPjwYYlVhkjjUPARIgMDBgzAvXv38Omnn2LkyJGYPn06cnJy2C5LZSl6i+/333/H8OHDsXXrVixatIgmsTQTBR8hMqKuro4ZM2YgJSUFZmZm6N69O9asWYOKigq2S1MplZWVSExMhKurK9ulNJpIJEJQUBC+//57XLlyBWPHjmW7JKVAwUeIjOnr62PdunW4e/cuYmNj4eDggJCQEBr/ayHx8fHo3Lkz9PT02C6lUcrLyzF58mT89ddfuHv3Lnr06MF2SUqDgo+QFmJra4vQ0FAcPHgQ69evh6enJ6Kjo9kuS+kp4oXrWVlZ8PLygpaWFsLCwugSGSmj4COkhXl7eyM6OhrTp0/HBx98gGnTpiErK4vtspSWol24Hh0djd69e2PChAk4ePAgTWKRAQo+Qligrq6O6dOnIyUlBZaWlnB2dkZwcDDKy8vZLk3pKFKL78SJExg5ciS2b9+OwMBAmsQiI7RyCyFyICMjAwEBAbhz5w7Wr1+PyZMn04eeFLx+/Rrt27dHYWEh1NXV2S6nTm8nsRw6dAhnz56l7a9kjFp8hMgBa2trhISE4NixY9i0aRM8PDwQFRXFdlkKLzo6Gm5ubnIdemVlZZg0aRKuXbuGqKgoCr0WQMFHiBzp378/oqKiMHv2bIwfPx5Tp07F8+fP2S5LYcn7+F5mZiY8PT2hq6tLk1haEAUfIXJGTU0Nn376KZKTk2FrawsXFxcEBQWhrKyM7dIUjjyP7729qH7y5Mk4cOAA7e3Ygij4CJFTenp6WLVqFWJjY5Gamgp7e3scOXIEIpGI7dIUAsMwctviO3bsGPz8/LBr1y788MMPNJ7bwmhyCyEK4u+//8aCBQsAAJs3b0bfvn1Zrki+ZWRkwMPDA1lZWXITLCKRCMuXL8fRo0dx9uxZdO/ene2SVBK1+AhREB4eHrhz5w7mzZuHSZMm4aOPPsLTp0/ZLktuvW3tyUvolZWVYcKECQgPD0dUVBSFHoso+AhRIGpqapgyZQqSk5Nhb28PV1dXLFu2DKWlpWyXJnfkaXzv+fPn6N+/PwwMDHDt2jWYmZmxXZJKo+AjRAHp6upixYoViI+PR0ZGBuzt7XHw4EEa/3uHvIzvRUZGok+fPpg6FKv9yQAACopJREFUdSr2799Pk1jkAI3xEaIEoqKi4O/vj+rqamzevBn9+/dnuyRWVVdXw8jICNnZ2dDX12etjqNHj2LBggXYv38/Ro0axVodRBKX7QIIIc3Xu3dv/P333wgJCcGUKVPQu3dvrF+/HjY2NmyXxooHDx6gQ4cOrIWeSCTC0qVLceLECYSFhaFbt26s1EFqR12dhCgJDoeDyZMnIykpCU5OTujZsycWL16MkpIStktrcWzuuF5aWorx48fj1q1biIqKotCTQxR8hCgZHR0dLF26FPHx8cjOzoadnR32798PoVDIdmkthq0d1589e4b+/fvD2NgYV69ehampaYvXQN6Pgo8QJWVhYYEDBw7gzJkz2L9/P3r27ImIiAi2y2oRbExsuXPnDvr06YNPP/0Ue/fuhYaGRos+P2k4mtxCiApgGAZ//PEHfvjhB7i5uWHDhg2wtbVluyyZKC4uhrm5OQoLC8Hj8VrkOQ8fPoxvv/0WBw4cwIgRI1rkOUnTUYuPEBXA4XAwadIkJCUlwc3NDb169UJAQACKi4vZLk3qYmJi0KNHjxYJPZFIhMDAQAQFBeH69esUegqCgo8QFaKtrY3FixcjMTER+fn5sLOzw549e5Rq/K+lLlwvKSnB2LFjcefOHURFRcHR0VHmz0mkg4KPEBVkbm6Offv24fz58zh8+DBcXV0RFhbGdllS0RLje0+fPkW/fv1gZmaGK1euwMTERKbPR6SLgo8QFebq6oqIiAgsX74cM2bMwJgxY/D48WO2y2oWWbf4bt++jT59+mD69On49ddfaRKLAqLgI0TFcTgcjB8/Ho8ePULfvn3Rt29ffPfddygsLGS7tEbLzMxEVVUVrK2tZXL+gwcPYuzYsdi/fz/8/f3lZgFs0jgUfIQQAICWlhYCAgLw4MEDFBUVwd7eHrt27YJAIGC7tAZ729qTdiAJhUL88MMPCA4ORkREBIYPHy7V85OWRcFHCJHQtm1b7NmzB5cuXUJISAhcXFxw5coVtstqEFlcuF5SUoIxY8YgOjoaUVFRcHBwkOr5Scuj4COE1KpHjx4ICwtDcHAwZs+ejdGjRyM1NZXtsuol7aXK3m5m265dO1y+fBmtW7eW2rkJeyj4CCF14nA4GDNmDB4+fAgvLy94eHhgwYIFeP36Ndul1SAUCnHv3j307NlTKue7desW+vbti5kzZ2LXrl00iUWJUPARQt5LU1MT3333HR49eoSKigrY29tj+/btcjX+l5SUhLZt28LY2LjZ5/rtt98wbtw4HDhwAPPmzaNJLEqGgo8Q0mBmZmbYtWsXrly5glOnTsHZ2RmXL19muywA0hnfEwqF+O6777BmzRrcuHEDQ4cOlVJ1RJ7QfnyEkEZzcnLClStXcO7cOXzzzTfo0qULNm7cCHt7e9Zqau74XnFxMT7++GOUl5cjKipKKi1HIp+oxUcIaRIOh4PRo0fj4cOHGDRoEDw9PTFv3jwUFBSwUk9zWnxPnjyBh4cHrKyscPnyZQo9JUfBRwhpFg0NDSxcuBCPHj2CUCiEg4MDtm7diurq6haroaysDI8fP4azs3OjH3vjxg14eHjgq6++wo4dO1psRwfCHtqWiBAiVQ8ePMDChQvx7NkzbNq0CcOHD5fZ5JD09HTs2bMHenp6+OOPP3D//v1GPX7fvn1YtGgRjh49isGDB8ukRiJ/KPgIIVLHMAwuXLiAhQsXwsbGBhs3bpTJ7gV37tyBp6cnuFwuBAIBDAwM8Ouvv2L8+PH1Pk4oFOL777/Hf//7X5w7dw52dnZSr43IL+rqJIRIHYfDwciRI/HgwQMMHz4cPj4+mDNnDvLz86X6PO7u7uDxeKisrIRQKASfz0e3bt3qfUxRURH8/PyQkJCAyMhICj0VRMFHyP+1dz+hUVwBHMd/O7ubnU10mxqzrCZSA6FZMVbQQw1CTIKICipYe2ktvQVMDWIS0NZSZMH0UDGUEEsQb2m1kFMFA8nBLIHSClra0maVSC2Nf9qYRWJq9l82vTSh6yZarVlnO9/P8c2bx5vL/ObNzHsPi8btduvQoUOKRCIyDENr1qxRR0eHEonEc2t/dsK6aZo6d+5cVpD986XWjRs3VFNTo4qKCvX19fETi00RfAAWXUlJiTo7OxUOhzUwMKDq6mpduHBBz+NLS319vSSppaVFu3fvzjgWjUYVCAQUDoc1ODiozZs36+DBg+rq6uInFhvjGx+AnOvr61Nra6vKysp06tQprVu37pnbunz5sg4fPqyhoSEZRuazfHt7u44fPy6n06mioiKdP39eW7du/a/dR54j+AC8EMlkUt3d3QqFQtq7d69CoZD8fv8Tz7s3GVfvlVFF7k5oIpaSz3QpGPDpzY3lKlniyWg/EAgoGo1KklauXKnh4WH5fL5FuybkB4IPwAsVjUYVCoXU09Ojo0ePqrm5WR6PJ6ve97/dV9fgiMLXxyRJ8VR67pjpMjQjqa6qVE1bKrV+VbHOnj2rxsZGpdNpmaapRCKh9vZ2HTlyJFeXBosi+ABYQiQSUVtbmyKRiE6ePKk9e/bMzf/r+eamTlyMKJaa1uPuWA6HZLqcOrYzqI/eatDt27fV0NCgHTt2qLa2VtXV1XI6nTm6IlgVwQfAUvr7+9XS0iK/36+Ojg79OPWSTlwc1lQy/eST/+Z1G3p/e1Dv1KxmZwVkIfgAWE4qldKZM2f0cfcXKtz1gWKpzNCbnnqg8YufKnbzOxlen17e8q6K1tZl1PG6nfqycZNeKy/OYc+RD5jOAMByXC6XDhw4oO2tHYpPZ4/0ov2fyeF0q7y5R8t3tWm8/7QSY79m1ImlpnV6cCRXXUYeIfgAWNK9ybjC18eyvumlEzE9vPa1imv3yyjwyly1VoWVr+vPny5l1JuZkS5dG9P4ZDyHvUY+IPgAWFLvldF5y1PRW3IYhtzLyubK3P4KJR8Z8UmSQ1Lv1fnbgX0RfAAsKXJ3ImPKwqx0ckoOT2FGmeEpVDoxlVU3lkorcufBovUR+YngA2BJE7HUvOWG26uZeGbIzcQfyijwLtBO7vYFRH4g+ABYks90zVvuWlammfS0ktFbc2WJP36Ru/SVBdphTU5kIvgAWFIw4JPHlX2LMgpMFVbV6P7Q50onYoqN/qyHI9+qaG19Vl3TZSi4Ymkuuos8QvABsKR9G8sXPLZsW5NmUgmNdr6te199opJtTSqYZ8Q3I2nfhoXbgT3N/y4BAF6w5Us82vJqqQaGf8+a0uD0LpX/jQ8fe77DIdVXlWYsXA1IjPgAWNh7dZUyXc+2tqbpcqqprvI59wj/BwQfAMtav6pYx3YG5XU/3a3K6zZ0bGeQ5cowL151ArC0/ZtWS9JT784wex7wKBapBpAXfhi9r9ODI7p0bUwOKWPh6tn9+OqrStVUV8lID49F8AHIK+OTcfVeHVXkzgNNxJLymW4FVyzVvg3l/MiCf4XgAwDYCj+3AABsheADANgKwQcAsBWCDwBgKwQfAMBWCD4AgK0QfAAAWyH4AAC2QvABAGyF4AMA2ArBBwCwFYIPAGArBB8AwFYIPgCArRB8AABbIfgAALZC8AEAbIXgAwDYCsEHALAVgg8AYCsEHwDAVv4CGjhAbsvsjZYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nx.draw(graph, with_labels=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we aggregate the all the transitions into *incoming traffic* and *outgoing traffic*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "incoming traffic: [ 266.    0.  247.  385.  155.  583.  217. 1147.]\n",
      "outgoing traffic: [376. 380. 386. 365. 370. 385. 376. 362.]\n"
     ]
    }
   ],
   "source": [
    "traffic_in = transitions.sum(axis=0)\n",
    "traffic_out = transitions.sum(axis=1)\n",
    "\n",
    "print(\"incoming traffic:\", traffic_in)\n",
    "print(\"outgoing traffic:\", traffic_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Estimating transitions using ChoiceRank\n",
    "\n",
    "ChoiceRank can be used to recover the transitions on the network based only on:\n",
    "\n",
    "1. information about the structure of the network, and\n",
    "2. the (marginal) incoming and outgoing traffic at each node.\n",
    "\n",
    "ChoiceRank works under the assumption that each node has a latent \"preference\" score, and that transitions follow Luce's choice model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = choix.choicerank(graph, traffic_in, traffic_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can attempt to reconstruct the transition matrix using the marginal traffic data and the parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "est = np.zeros((n_items, n_items))\n",
    "for src in range(n_items):\n",
    "    neighbors = list(graph.successors(src))\n",
    "    if len(neighbors) == 0:\n",
    "        continue\n",
    "    probs = choix.probabilities(neighbors, params)\n",
    "    est[src,neighbors] = traffic_out[src] * probs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True transition matrix:\n",
      "[[  0.   0.   0.   0.   0.   0.   0. 376.]\n",
      " [  0.   0.   0.   0.   0. 380.   0.   0.]\n",
      " [  0.   0.   0.   0.   0.   0.   0. 386.]\n",
      " [148.   0.   0.   0.   0.   0. 217.   0.]\n",
      " [118.   0.  47. 135.   0.  70.   0.   0.]\n",
      " [  0.   0.   0.   0.   0.   0.   0. 385.]\n",
      " [  0.   0.  88.   0. 155. 133.   0.   0.]\n",
      " [  0.   0. 112. 250.   0.   0.   0.   0.]]\n",
      "\n",
      "Estimated transition matrix:\n",
      "[[  0.   0.   0.   0.   0.   0.   0. 376.]\n",
      " [  0.   0.   0.   0.   0. 380.   0.   0.]\n",
      " [  0.   0.   0.   0.   0.   0.   0. 386.]\n",
      " [149.   0.   0.   0.   0.   0. 216.   0.]\n",
      " [117.   0.  52. 127.   0.  74.   0.   0.]\n",
      " [  0.   0.   0.   0.   0.   0.   0. 385.]\n",
      " [  0.   0.  91.   0. 155. 130.   0.   0.]\n",
      " [  0.   0. 105. 257.   0.   0.   0.   0.]]\n",
      "\n",
      "Difference:\n",
      "[[ 0.  0.  0.  0.  0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.]\n",
      " [-1.  0.  0.  0.  0.  0.  1.  0.]\n",
      " [ 1.  0. -5.  8.  0. -4.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.  0.  0.]\n",
      " [ 0.  0. -3.  0. -0.  3.  0.  0.]\n",
      " [ 0.  0.  7. -7.  0.  0.  0.  0.]]\n"
     ]
    }
   ],
   "source": [
    "print(\"True transition matrix:\")\n",
    "print(transitions)\n",
    "\n",
    "print(\"\\nEstimated transition matrix:\")\n",
    "print(np.round_(est))\n",
    "\n",
    "print(\"\\nDifference:\")\n",
    "print(np.round_(transitions - est))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
